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Abstract. The paper presents a globally convergent algorithm for solving coefficient inverse 
problems. Being rooted in the globally convergent numerical method (SIAM J. Sci. Comput., 
31, No.l (2008), pp. 478-509) for solving multidimensional coefficient inverse problems, it has two 
distinctive features: the new iterative and refinement procedures. These novelties enhance, some- 
times significantly, both the spatial and contrast resolutions. The computational effectiveness of 
the proposed technique is demonstrated in numerical experiments with two applied coefficient in- 
verse problems: electromagnetic or acoustic frequency sounding and electrical prospecting of layered 
media. The Slichter-Langer-Tikhonov formulation is exploited as a mathematical model of the latter. 
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1. Introduction. The paper addresses a challenging problem of constructing 
globally convergent algorithms for Coefficient Inverse Problems (CIPs), i.e., the in- 
verse problems for the partial differential equations where one needs to determine one 
or several coefficients from a given solution (or its functionals) on some manifolds 
(see, e.g., [3J, [7J, [TS] and references cited there). In this paper, the global conver- 
gence is understood in the sense that a reconstruction algorithm produces a sequence 
of iterates that converge to the sought coefficient or to its "good" estimate from an 
initial approximation, which is not necessarily close to this coefficient. 

Due to both the nonlinearity and ill-posedness of CIPs, providing the global 
convergence is important for many, if not for all, reconstruction algorithms. Being 
applied to a CIP, the nonlinear least squares approach normally results in a nonconvex 
minimization problem. Clearly, under such conditions, the traditional numerical tech- 
niques, such as the gradient or Newton-like methods, may be not applicable or they 
may fail to converge to the true solution. On the other hand, the methods of global 
optimization are extremely time consuming and heuristic. Thus, the development of 
globally convergent algorithms is of particular interest to quantitative imaging in de- 
fence science, geophysics, medical diagnostics, non-destructive testing and evaluation, 
etc. 

The goal of this paper is twofold. First, this is to present a globally convergent 
algorithm for solving some CIPs in one dimension. Second, this is to demonstrate 
its computational effectiveness by applying it to frequency sounding and resistivity 
prospecting of layered media. Introduced by Tikhonov |26j and Cagniard [4j in explo- 
ration geophysics, the method of frequency sounding became one of the most popular 
techniques used in a variety of applications (see, e.g., survey [9]). The term "resistivity 
prospecting" is associated with geophysical prospecting techniques utilizing measure- 
ments of the voltage potential produced by currents injected from the Earth's surface 
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into the ground (see, e.g. [30] ). Unlike the Calderon formulation which is widely 
used in electrical impedance tomography with multiple boundary measurements, we 
utilize the Slichter-Langer-Tikhonov mathematical model [T3] , [3D] , [37] with a single 
boundary measurement. Based, in general, on three dimensional mathematical mod- 
els, these CIPs allow, nevertheless, for one dimensional formulations making, by the 
same token, more descriptive the underlying ideas of the proposed algorithm. On the 
other hand, one dimensional quantitative imaging is widely used in many applications, 
such as studying Earth interior [3T] and seafioor resources [5], as well as in studying 
macromolecular biological interactions [22j . nuclear magnetic resonance ID imaging 
[16) . mine countermeasures [121 . etc. 

The globally convergent numerical method for solving the hyperbolic CIPs was 
proposed and developed in [2J, [3J, [7], [ID], [ZD, H2- ^ n the authors utilized the 
Carleman estimates and constructed some strictly convex least squares functionals 
at each step of a layer stripping procedure with respect to the spatial variable. In 
H], 0, [ID], [H], [U] the layer stripping procedure with respect to the real- valued 
parameter of the Laplace transform combined with an iterative procedure for treat- 
ment of the so-called tails was developed and applied to solving jiD (n = 1,2,3) 
CIPs. Being rooted in the globally convergent numerical method [2], the proposed 
algorithm possesses two new features. First, in the iterative procedure for solving an 
overdetcrmined boundary value problem for the nonlinear second-order differential 
equation the successive approximations are constructed for every parameter from the 
Laplace transform. Second, we refine the iterative solution by matching the inte- 
rior field generated by this solution with the solution of the forward problem in the 
region of interest unlike the traditional refinement by matching the boundary data. 
This technique differs significantly from a refinement procedure (see [3J) based on the 
adaptive FM/FD methods [T]. 

The paper is formatted as follows. In the next section, a generic coefficient in- 
verse problem is formulated to exemplify our approach. In section [3] we describe 
the proposed algorithm. In sections |4] and [5] we state convergence results for both 
the successive approximations produced by the proposed algorithm and their refine- 
ments. In section[6]we conduct the numerical study to demonstrate the computational 
effectiveness of the proposed algorithm. We conclude our investigation in section [Jj 

2. Problem formulation. We introduce a generic coefficient inverse problem 
as follows. Let C i? 3 be a convex bounded domain with a sufficiently smooth 
boundary d£l. Let the function /i(x) > 1 be such that /i g C 2 (R 3 ), fi(x) = 1 for 
x £ i? 3 \fi. Consider the Cauchy problem for the hyperbolic equation 



In particular, if the problem (|2.1[) - (|2.2j) models the propagation and scattering of 
the electromagnetic or acoustic field, then = c~ 2 (cc), where c{x) is the speed of 
electromagnetic or acoustic waves. We emphasize that we do not impose the smallness 
assumption on the coefficient /i. 

Generic Coefficient Inverse Problem (GCIP). Let the function U(x,t) be 
the solution of the Cauchy problem (2.1]) - I12.2]) . Given a single measurement U(x,t) = 
f{x 1 1) on the boundary dfl for a fixed source position xo ^ f2 and for all t G (0, oo), 
find the coefficient jj,{x) in tt. 

To our knowledge, there is no global uniqueness result for the GCIP available in 
the mathematics literature. We assume below uniqueness of the solution of GCIP, 



(2.1) 
(2.2) 



(i(x)Utt - AU = in R 3 x (0, oo), 
U{x,0) = 0, U t (x,0) =5(x-x ). 
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because we will utilize the generic inverse problem only for the purpose of describing 
the globally convergent algorithm. Consider the Laplace transform of the function 
U(x, t) for s > s, where s is a positive constant, 

/•OO 

u(x,s)= / U{x,t)e~ st dt, s > s = const. > 0. 
Jo 

It follows from [3] (see Theorem 2.7.1) that for all s > s, where s is sufficiently large, 
the function u(x, s) is the unique solution of the problem 

(2.3) Au — s pb(x)u = —5(x — xq), x £ R 3 , 

(2.4) lim u(x,s) = 



satisfying the following conditions 



(2.5) u(x, s)>0yx^x o ;u (x, s) = GXp( ^ X f o|) + u (x, s) , u £ C 2+a (R 3 ) 

Att\x — xq\ 

for all a £ (0, 1). Furthermore, for any s > the theorem 2.7.2 in [3] ensures existence 
of a unique solution of the problem (|2 . 3|) , (|2.4I) satisfying conditions [23]) . In this paper, 
C k+a are Holder spaces, where k > is an integer and a £ (0, 1). 

For simplicity, we describe the proposed algorithm in one dimension. Its general- 
ization to two and three dimensions is straightforward. In one dimension, SI = (0, 1), 
and the source is located outside of this interval, i.e., xq £ {x < 0} . Also, we assume 
that the function fi(x) satisfies the following conditions 

(2.6) n £ C a (R) ,fi(x) £ [1, m],[i(x) = 1 for x £ (0, 1) , 

where the number m > 1 is given. Then a single measurement is given by {7(0, t) = 
fit), and the one dimensional analogue of the problem (|2.3[) - (|2.4[) acquires the form 

(2.7) u xx — s 2 /i(x)w = — S(x — Xq), X £ R, 

(2.8) lim u(x,s) = 0. 

It follows from [12] (see the theorem 3.1) that the Laplace transform of U (x, t) satisfies 
the problem (|2.7j) - (|2.8j) for all s > s for a sufficiently large s > 0, and for all s > there 
exists a unique solution u(x, s) of the problem (|2.7[) . (|2.8[) . It satisfies the analogue 
of (EJJ 

(2.9) u(x, s) > 0,Vx £ R; u (x, s) = u (x, s) +u(x,s) ,u ' G C 2+a (R) , 
where 

exp (— s \x — Xq\) 



Uq(x,s) 



2s 



is the solution of the problem (|2.7I) - (|2.8[) with ^(a;) = 1. Let ip{s) be the Laplace 
transform of fit). Since n(x) = 1 for x < 0, then given y(s), one can uniquely 
determine the function u x (0,s) = sip(s) — exp(sxo) := (p(s) (see, e.g., [T2"]). 
Now we introduce a new function 

(2.10) «(x,a) = -ln(-). 
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Then the problem (|2.7I) - (|2.8[) is reduced to the problem 

(2.11) v xx + s 2 vl - 2sv x = (i(x) - 1 in (0,1), 

(2.12) v(0,s) = v o (s), v x {0,s) = Vl (s), 



where 



In ip(s) - In 2s x , , 2 e 
vo{s) = n 1 > v n s ) = — 



s ' s s 2 ip(s) 

Introducing one more function q[x 1 s) — d s v(x,s), we eliminate the unknown cocffi- 
cient fj,(x) from equation (|2.11l) by differentiating (|2.1ip . (|2.12l) with respect to the 
s-variable in accordance with [10] and [12]. As a result, we obtain the problem 

(2.13) q xx + 2s 2 q x v x - 2sq x = 2v x - 2sv% in (0, 1), s G [s, s], 

(2.14) q(0,s)=v' (s), q x (0 ) s)=v' 1 (s), 

(2.15) g x (l,s) = 0. 

Note that the condition (|2.15l) follows directly from (|2.8p . because fi(x) = 1 outside 
of the interval (0, 1) and 



(2.16) v(x, s) = — j q (x,t) dr + v (x, s) . 

where 

v(x,s) = — I q(x,T)dr 



3. Description of the algorithm. It follows from the previous section that if 
both functions q(x, s) and v(x, s) are known, then the solution of GCIP can be found 
in the closed form from the equation (|2.11[) . However, these functions are unknown. 
If one knows an approximate solution (q, v) to this problem, then the approximate 
solution fl(x) of the generic inverse problem is given by 

(3.1) p,(x) fa v xx + s 2 v 2 x - 2si x + 1. 

Thus, determining the approximate solution of the problem (|2.13|) - ()2.15j) plays the 
crucial role in developing the reconstruction algorithm. 

Note that if the function v(x, s) is given, then the problem (|2.13l) - (|2.15l) is linear 
with respect to q(x, s). It may seem, therefore, that in order to produce an approx- 
imate solution to this problem, one may solve a traditional boundary value problem 
for the second-order ODE (|2.13ll with two boundary conditions (|2.14p and (|2.15|) . 
However, since q(x, s) depends on v(x, s), one needs to iterate both functions q and 
v . In this case, the sucessive approximations may not converge to the appropriate 
approximations of q(x, s) and v(x, s). This fact was observed in |10] as well as in our 
numerical experiments when attempting to solve the traditional BVP. Therefore, the 
method of quasi-reversibility (see, e.g., [14 ), which is well suited to handle overde- 
termined boundary value problems, has been used for the approximate solution of 
the problem (|2.13[) - (I2.15I) . Let the function v(x,s) be fixed for any fixed parameter 
s € [s, s), and the left- and right-hand sides of the equation (|2.13p be denoted as 
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L(q)(x) and F(x). Then the method of quasi-reversibility consists of minimizing the 
Tikhonov functional 

(3-2) T x (q) = \\L(q)(x) - F(x)\\l M + X\\qf H2{0>1) , 

where A > is the regularization parameter, subject to the boundary conditions 
(|2. 14|) - (|2.15p . It was shown in [12] that this functional has a unique minimizer. 
Furthermore, since this functional is the sum of squares of two expressions, which 
depend linearly on the function q , it is strictly convex. This means that any numerical 
method, such as the gradient or Newton-like one, allows for determining numerically 
that unique minimizer q\ of the functional T\(q) for every s G [s, a]. However, it 
should be emphasized that in accordance with the theory of ill-posed problems (see, 
e.g., )29] ), the regularization parameter A must be chosen consistently with the level 
of noise 6 in the "measured" data (p(s), i.e., such that \\ip — <f>\\c[s,s] < S, in order 
to ensure the regularizing property. A priori parameter choice A(<5) = 5 2v , where 
v G (0, 1) (see, e.g., (25] [29]) is recommended for use. 

The regularization procedure described above allows us to solve numerically the 
problems (|2.13|) - (|2.15|) only if the function v(x,s) is known. However, in reality 
both functions q(x, s) and v(x, s) are unknown. To overcome this difficulty, we deter- 
mine the successive approximations of these functions as follows. We first utilize the 
asymptotic behavior of v(x, s). It follows from [10] (see Lemma 2.1) that there exists 
a function p(x) G C 2+Q [0, 1], such that 

v (x, s)^^- + (J^j , s ^ oo , x G [0, 1] . 

Since q{x, s) — d s v(x, s), we obtain for sufficiently large s 

(3.3) v(x, s) fa ,q(x, s) pa 

s s 

Because of asymptotics (|3.3p , we choose s = s and substitute these representations in 
(|2. 13|) - (|2.15p . As a result, we obtain the following problem 

(3.4) p xx (x) = in (0,1), 

(3.5) p(0) = ~s 2 v' (s), p' (0) = -s 2 Wl (s), p'(l) = 0. 
This problem is also overdetermined. 

Initialization. We solve numerically the problem (|3.4p ~ (|3.5l) by the method of quasi- 
reversibility. Denote pqrm (%) its solution. Then we determine the initial approxima- 
tion of the function v(x, s) as 



(3.6) 



vW(x,s) = Pqrm{x) , xe[0,ll, «gM- 



Clearly, the initial approximations of the functions u(x, s) and q(x, s) will be uq(x, s) 
and q°(x,s) = -Pqrm(x)/s 2 . 

Step 1. Beginning with v (x, s), we start the iterative process. Assume that 
the A;th successive approximation v^(x, s) has been determined for all parameters 
s G [s, s). Then, given v^(x : s), we solve the overdetermined problem for 

d 2 q {k+1) (x, s) + 2s 2 d x q {k+1) (x, s)d x v^ (x, s) - 2sd x q {k+1) (x, s) 
= 2d x v( k \x,s)-2s(d x v( k \x,s)) 2 in (0,1), 

q {k+1) (0,s) = tp' (s), d x q( k+1 \0, S ) = ^i( S ),^ g ( fc+1 )(l, S ) = 
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for each fixed parameter s € [s,s) and find the successive approximation q k+1 (x, s). 
Step 2. We use the formula (|2.16p to update the function v(x, s) for each parameter 
s £ [s, s) as follows. 



(3.7) v {k+1) {x,s) = -J q {k+1 \x,v)dv + v (k \x,s), 

(3.8) d x vi k+1) (*, s ) = -[ d -1 ik+1) v ) dv + d - v(k) *)> 

(3.9) 9> (fc+1) (x,s„_i) = - / d 2 x q {k+1 \x,v)dv + dlv {k \x,s). 
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Step 3. Let s* € [s, s) be a certain number which is independent on k. We update 
the coefficient as 

(3.10^+1(3:) - dy k+1 \x,s*) + (s*) 2 (d x v^ k+1 \x, S *)) 2 - 2s*d x v^ k+1 \x,s*) + 1, 
where x € (0, 1), and set (Xk+i{%) — 1 for x — and x = 1. 

iSiep ^. Update the function u(a:, s) and its derivatives by solving the forward problem 

(3.11) d^ k+1 \x,s)-s 2 fi k+1 (x)u(x,s) = -6(x-x ), ie(fl,l), xoi (0,1), 

(3.12) lim u (A;+1) (a;,s) = 

| a; I — >oo 

and determining the functions 

(3.13) v ^\xrs) = ^^il) _ £^ - ^1 

s s s 

(3 14) 9 (x s) d * u(k+1) ( x >s) I 

for x e [0, 1]. It follows from (l3Jl) - ([3~9)) that knowledge of the function u( fc+1 ) [x, s) 
and its derivatives is crucial in reconstructing n(x) via the formula (|3.10l) . 
Stopping criterion. The iterative process continues until a stopping criterion is ful- 
filled. If the level of noise 6 in the data ip(s) is given or it can be estimated, then it 
is natural to stop the process at the iterate 

(3.16) k stop = min{/c : \\d x Un k (0,s) - <£(s)||l 2 [s,s] < ^} , 

where the function it Mfc [x, s) is determined by solving the forward problem (|2.7j) - (|2.8[) 
with fj,(x) — fj>k(%) for every successive approximation fj-k(x). Note that since the 
norm in (|3.16[) represents the residual functional on Hk{x), the number k stop plays 
the role of the regularization parameter. In this case, the closeness of [ik Etov {%) to the 
true coefficient n(x) needs to be estimated. Also, we note that in computing we deal 
with the discrete functions defined on grids with respect to both variables x and s. 
In this case, the derivatives are approximated by their finite difference analogous and 
the integrals in (|3.7I) - (|3.9[) - by appropriate quadratures. 
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4. Convergence of successive approximations. Let s > 0, s > 1. Let ^*(x) 
be the unique solution of the GCIP satisfying conditions (I2.6P and corresponding to 
the data y(s), s € [s, s]. In the convergence analysis we assume that all successive 
approximations /ifc(x) obtained by the algorithm described in the section [3] also 
satisfy the conditions (|2.6I) . In addition, following (|2 . 10[) and (|3.3[) . we assume that 
the function v (x, s) = s~ 2 In (u/uq) can be represented as 

(4.1) u(o; ) s) = ^,Vs>s, 

s 

where the function p 6 C 2+a [0, 1]. Let /i > be the largest distance between two 
adjacent nodes in a grid with respect to the s-variable. Let x & (0, 1) be a given real 
number, such that \\<p — <p\\c[s,s\ < X, where <p(s) is the single boundary measurement. 
Denote r\ = max(<5, h). Then the following analogue of Theorem 6.1 in [12] takes place. 

Theorem 4.1. If the initial approximation v^(x,s) in the proposed algorithm 
is determined as in 13. 6\) . and the regularization parameter A in the Tikhonov's func- 
tional \3. £j) is chosen as A = % 2 , then there exists a constant C = C(xo,to,s), such 
that for r\ chosen from 

1 

n < 



the following inequality is fulfilled for all k £ [1, k stop \ 

(4-2) IMz)-M*(*)lk(0,i) <if, 7€ (0,1), 

where the number"/ does not depend on k,h,x, fJ-k(x) and /z*(x). 

Proof. Since the arguments are analogous to those used for the proof of Theorem 
6.1 of [12], we only outline the scheme of the proof. First, using p.4[) - (|3.6|) . (|4.1[) and 
Lemma 6.1 in [12], it is shown that 



(4.3) max (x, s) — v (x, s) 

se[s,s] 



< Of]. 

H 2 (0,1) 



Recall that to is the constant in (|2.6p . Consider the solution ujt(x, s) of the problem 
([2~7]) . PT5|) for the case /i (x) = to, 

exp(-s % /m|o; - x |) 
2s 

Then it follows from [12] 

(4.4) u m (x, s) < u(x,s) <uq (x, s) , Vx € [0, 1] , Vs e [s, s] . 



Using (|4.4I) . as well as estimates of a minimizer of the functional (|3.2p indicated in 
Lemma 5.2 in [T2], one can estimate norms g« — g n , n € [1,-AH. Next, 

fl" 2 (0,l) 



using (|3.7p - (|3.15p and (|4.3[) . one can obtain the estimate (|4.2p for k = 1. Next, 
using the latter estimate as well as (14.41) . one can obtain an analogue of (|4.3|) for 
max sS [s s] ||u < - 1 ' (x, s) — v (x, s)|| ff2 ( jy Finally, we can continue this process in a sim- 
ilar way until it stops at k = k s t op . □ 

The main consequence of this theorem that it guarantees a good approximation 
of the true solution pL* regardless a priori knowledge of a small neigborhood of this 
solution. 



8 



M. V. KLIBANOV, AND A. TIMONOV 



5. Refinement of //fc sto . The theorem 14.11 guarantees an appropriate closeness 
of successive approximations to the true coefficient only if numbers x an d h are suffi- 
ciently small. In practice, however, it is more typical than exceptional that the level of 
noise \ i n the experimental data is fixed and it is not small enough. In this case, the 
approximation Hk, top (x) may not possess high spatial and contrast resolutions, i.e., 
the quantitative imaging may be problematic. To alleviate this restriction, we propose 
to take advantage of knowledge of not only ^k stop {x) but also the field v^- kstop \x, s) 
everywhere in [0, 1]. Indeed, we observe that given i)( fc <"°*>) (x, s), one can determine 
from (|2.10[) the function 

(5.1) w (fcstoJ, ^(x, s) = u (x, s) exp [s 2 v^ katop ^ (x, s)}. 

Choose a parameter s* G [s, s) and denote u(x, s*) = u^ kstopS> (x, s*) the interior data. 
We observe that although the function u(x, s») may be close to the solution of (|2.7I 
)- (|2.8p . it does not satisfy this problem. This observation makes it possible matching 
the interior field u(x, s*) with the function u M (x, s») with the hope of improving the 
approximation fik sto {%)■ Here, u M (x, s),s £ [s,s] is the solution of the problem (|2.7p ~ 
(|2.8|) that corresponds to the function fi(x) satisfying conditions (|2.6j) . 

Although the results established below remain valid in the Banach spaces, we 
consider the finite-dimensional Hilbert spaces in view of their use in scientific com- 
puting. Let H be such a space. As an example, we indicate the subspace of the space 
Li (R) spanned by a finite number of piecewise linear finite elements vanished out- 
side of the interval (0,1). For every function \i (x) satisfying conditions (|2.6[) denote 
H{x) — 1 = ~p{x) £ L2 (-R). Let G C H be the set of such functions. Clearly, G is 
bounded in the space H . Consider the map F : G — > L2 (0, 1), where 

F (71) = u M (x, s*) ,x 6 (0, 1) ,V~p(x) = li{x) — 1 € G. 

It follows from (|2.7[) that the operator F is one to one. Denote K = F (G) the range 
of F. Since G is a compact set in H, then the Tikhonov's theorem [3], [29] implies 
that the inverse operator F^ 1 : K —> G is continuous. Let luf (z) , z £ {z > 0} be the 
modulo of continuity of the operator on the set K. Then 

(5.2) llMx -M2ll.fr < (ll F (Mi) - f (M2)IIl 2 ( ,i)) y~Pn~p2 e G. 

In the GCIP both functions fi* (x) and it*»(a;,s*) are supposed to be unknown. In- 
stead, their approximations fi = ^k stop {x) G G (see (|4.2I) ) and u{x, s*) £ L2 (0, 1) are 
given, such that 

(5.3) (x, s*)-u (x, s*) ||i 2 ( ,i) < Q, 

where g = Crf amraa , C = const > 0. We are interested in finding a new approxima- 
tion of n*(x) that would be more close to it than the approximation fik 3to ( x )- 
Consider the Tikhonov functional 

(5.4) T X (7Z) = \\F (/I) - fi (x, s,) ||| 2(0)1) + A||m - JIk 3top \\h, 

where ~p, J^k Btop € G. Since G is a bounded compact set in the finite dimensional space 
H, then it follows from the Weierstrass theorem that for any A > there exists a 
minimizer j-i\(x) £ G of the functional (|5.4I) . i.e., 

fi\(x) = argmin {T\(u) :Jl£G) . 
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Moreover, one can show that the map F has the Lipshitz-continuous Frechet deriva- 
tive. Then, it follows from Theorem 1.9.1.2 in [3J that the Tikhonov's functional T\(jl) 
is strictly convex in a sufficiently small neighborhood of the function jf . This property 
allows for establishing convergence of (x) to \x* (x) . For any positive real number a 
and any x G H we denote U a (x) = {z G H : \\x — z\\ < a} the a - neighborhood of the 
element x. Then the following theorem takes place. 

Theorem 5.1. Let A = X(g) = g 2u ,v G (0,1/2), and £ G (0,1) be an arbitrary 
number. Then there exists a sufficiently small number go = go (H, m, v, £) such that 
for all g G (0, go) 



(5.5) 



Ma(<5) - I 1 



< 



H 



If v G (0,1/4), then there exists a sufficiently small number ~g = ~g(H,m,is) G (0,1) 
such that for every g G (0, <5) the functional j5.4\ ) is strictly convex in the neighborhood 
U g 3v (Jl*). If the solution /ifc stop is so accurate that — ~pk stop \\h < g 3l, /3, then 
~fi\( e ) € U g 3v / 3 (/I*), and any gradient or Newton-like method of minimization ofT\(jl) 
with the initial approximation ^k stop converges to A*A(g) ■ 

Proof. Since the arguments established in [5] and in theorems 1.8 and 1.9.1.2 of 
[3J work for this proof, we only prove (|5.5[) . We have 



-FQ?) 



< 



£2(0,1) 

Hence, (|5.3p implies that 



f((J.x(S)) ~u{x,s*) +\\F{fi*)-u{x,s^)\\ L2(01) 

* ' -^2 (Oil) 



(5.6) 



(iix(s))-F(r) 



< 



£2(0,1) 



F 



(ma(5)J -«(x,s*) 



£2(0,1) 



We now estimate the first term in the right hand side of (|5.6p . Since g is sufficiently 
small, g 2 < g 2v . Hence, we have 



F 



(^A( e )) ""( x ' s *)|| i2(01) < T A(e)(MA( e )) < r A(5)(M*) 



II/'" 



Then we obtain 



(MA (e) ) - l} < u ^+\\r--pk st j 2 H+Q < w^fi+ii 



Using (|5.2p . we obtain 



(5.7) 



Mau 



h < u F (2g^l + \\r-Jl kst J 2 H ) 



We first consider the option Hk stop A 1 *- Since lim^ + u f (z) = 0, then we can find 
such number go — go {H, m, v, £) G (0, 1) that for all g G (0, go) 



(2g^i + \\r -ji ksi j 2 H ) < eiir --p kat j H . 



(5.8) 



Combining (|5.7p and (|5.8p . we obtain the inequality (|5.5p . The second option can be 
considered by analogy. □ 

Note that although the theorem 15.11 does not allow for the quantitative estimate 
of the distance between fi* and /^a, it justifies the use of the Tikhonov regularization 
with the interior data for refining fik ato ■ 
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2 4 6 6 10 

Frequency 

Fig. 5.1. The simulated electromagnetic frequency dataTp(s) in the interval [1,10]. The noise- 
less data is shown by bullets and the data corrupted by noise at the level of 5% is shown by asterisks. 

6. Numerical study. In this section we perform some numerical experiments 
in order to demonstrate the computational effectiveness of the proposed algorithm. 
For this purpose, we show how it can be applied to frequency sounding and electrical 
prospecting of layered media. 
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Fig. 6.1. Reconstructions of the relative electrical permittivity of a single mid- contrast land 
mine in the air: (1) without refinement (dashed), (2) with refinement using the boundary data 
(asterisks), (3) with refinement using the interior data (bullets). The original distribution is shown 
by a solid line. 

6.1. Frequency sounding. Without loss of generality, we consider frequency 
sounding of layered media and exploit its mathematical model developed in [M] . In 
accordance with this model, a pulsed plane wave, which is normally incident at the 
plane z — 0, propagates through an inhomogeneous layer whose material property 
depends continuously on one variable z only and it is supposed to be constant outside 
of the interval (0, L). The propagation and scattering of such a pulsed plane wave 
can be described by the following initial boundary value problem 

(6.1) c- 2 (z)U tt -U zz = 0, z>0, t>0, 

(6.2) U(z ) 0) = U t (z,0) = O i 

(6.3) U(0,t) = S(t). 
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If the incident wave is electromagnetic and a medium is dielectric, the scalar quantity 
U(z, t) can be interpreted as an appropriate component of the electric field propagat- 
ing with the speed c{z) = \ j \J (1o(j,Eoe(z). Here, e{z) and fi are relative permittivity 
and permeability of the inhomogeneous layer and £o and A*o are absolute permittivity 
and permeability of vacuum. If the incident wave is acoustic, the quantity U(z, t) can 
be interpited as the acoustic potential and c(z) is the sound speed. 




Fig. 6.2. Recovering the relative electrical permittivity of a single mid- contrast land mine from 
the perturbed data (the level of noise is 5%) without refinement (asterisks) and with refinement 
(bullets). 

Inverse Problem of Frequency Sounding. Let the function U(z,t) be the so- 
lution of the initial boundary value problem I6.1\) - Wlfy . Given the function J7 Z (0, t) = 
*5f(t),t > and parameters cq > O.L > ; such that c{z) = Co for z > L, find the 
variable speed c(z) on (0,L). 

The traditional approach to this inverse problem consists of reducing the wave 
equation (|6.ip to the Hclmholtz equation by performing the Fourier transform and 
solving it by some globally convergent algorithms. To the author's knowledge, there 
are three groups of such algorithms. The algorithms of the first group are based on 
the so-called trace (asymptotic) formulae (see, e.g., [5]). The algorithms of the second 
group utilize the nonlinear Riesz transform allowing for reducing to an equivalent 
Volterra integral equation [23J . The algorithms of the third group are based on the 
convexification method [7] . Unlike the traditional approach, in [2] , [3J , [TU] , [H] , [H] , 
and [24) the Laplace transform was applied to an initial-boundary value problem for 
the wave equation. It was resulted in several other globally convergent algorithms. 

Following the latter approach, we apply the Laplace transform 

y* oo 

u(z, v) = / e~" t U(z,t)dt, z > 0, v > vq 
Jo 

to the problem (I6.ip - (|6.3p and introduce the dimensionless variables x = z/L,s = 
uL/co,n(x) = cq/c{Lx). Denoting u(x, s) — u(Lx,cos/L), we obtain the dimension- 
less two-point problem for the second-order ordinary differential equation 

(6.4) u xx (x, s) — s 2 n 2 (x)u(x, s) = 0, i£ (0,1), s G [s,s], 

(6.5) u(0,s) = l, 

(6.6) u x {l,s) + su(l,s) = 0. 
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It should be emphasized that if we consider the equation (|6.1|) for any z € (—00,00) 
together with the condition Ut(z, 0) = S(z — z ),z n < 0, then we obtain the GCIP 
formulated in the section [3J In this case, the transmission condition (|6.6[) follows 
from the radiation condition lim^^oo u(x, s) = and n(x) = 1 for x > 1. It follows 
from [24] that for a sufficiently smooth refraction coefficient n{x) the problem (|6.4I) - 
(|6.6p has a unique solution. Therefore, the inverse problem of frequency sounding can 
be reformulated as follows. 

Let the function u{x, s) be the solution of the problem jff.^p - jTTS]) that corresponds 
to n(x), but u x (0,s) is unknown. Instead, an approximation Tp(s),s € [s,s],s > of 
u x (0,s) and the real number S > are given, such that \\u x (0, s) — v(s)||i 2 [s,s] ^ 
Find an approximation of n(x) in (0,1). 




Fig. 6.3. Reconstructions of the relative electrical permittivity of a single low-contrast (left), 
mid-conrast (middle) and high- contrast (right) land mine in the air. The results of reconstruction 
are shown by asterisks (without refinement) or by bullets (with refinement). The original distribution 
is shown by a solid line. 

To generate the frequency sounding data p(s), we first simulated numerically the 
data for a given n(x). To provide a high accuracy of reconstruction, we represented 
the solution of the problem (I6.4p - (|6.6p in the form u = u e + u, where u e — e~ sx is its 
solution for n(x) = 1. Thus we solved numerically the boundary value problem 

(6.7) u xx (x, s) - s 2 n 2 (x)u(x, s) — -s 2 e~ sx (l - n 2 (x)), x 6(0,1), s € \s,s], 

(6.8) u(0,s) = 0, 

(6.9) u x (l,s)+su(l,s) =0. 

The second-order finite-difference analogue of this problem was numerically solved by 
a stable version of the elimination method (see [12]). Note that since the solution 
of this difference problem approximates the solution of the problem (|6.7[) - (|6.9p . the 
function <p(s) approximates u^O, s). The level of error of this approximation depends 
on several factors, such as smoothness of the refraction coefficient, grid, CPU, etc. In 
the numerical experiments, we utilized the uniform grids containing from 64 to 100 
nodes, and the refraction coefficient was assumed to be either smooth or piecewise- 
constant. Under such conditions, it was estimated from comparison with the model 
refraction coefficients that the level S varied approximately from 5 • 10~ 6 (if n(x) is 
smooth) to 5 • 10" 4 (if n(x) is piecewise-constant). 
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Fig. 6.4. Reconstructions of the relative electrical permittivity of the model clutter without 
(asterisks) and with (bullets) refinement. The original distribution is shown by a solid line. 



Since v = v — x/s and q = q + x/s 2 , where 
_ In (u + e~ sx ) x 



and q = d s v, 



the two-point problem for the function q has the form 

(6.10) q xx (x,s) + 2s(sv x - l)q x (x, s) = F(v x ; x, s), ie(0,l), s G \s,s], 

(6.11) g(Q,a) = 0, 

(6.12) q x (l,s) = 0, 

where F(v x ;x, s) — —2v x (sv x — 1). The boundary data is given by 

'Tp(s) 



(6.13) 



<l x (x,0) = tp(s) = d s 



In computing, the operator d s is understood in the sense of differentiating an interpo- 
lating cubic spline of Tp(s). We emphasize that the nonlinear problem (|6.10[) - f|6.13[) has 
two unknowns (v,q), and it looks similar to the problem (|2.13[) - (|2.15|) . Therefore, 
we applied the algorithm described in the section [3l Indeed, to start the iterative 
process, we let v^°'(x, s) — Pqum{ x )> w bere Vqrm i x ) ^ s the solution of the problem 



(6.14) 
(6.15) 
(6.16) 



Pxx = 0, xe (0,1), 
p(0)=p x (l) = 0, 
p x (0) = -s 2 ^(s) 



obtained by the method of quasi-reversibility. Suppose that after k iterates the ap- 
proximation y( k \x, s) is obtained. Then solving the problem (|6.10p - (|6.12[) with the 
boundary data ip( s ) by the quasi-reversibility method as described in the section [3l 
we obtained the corresponding approximation q( k \x,s). Given both these approxi- 
mations, we updated the refraction coefficient as follows. 



n 2 k+1 (x) = 1 + (G k (x) - 2s*H k (x) + s*H 2 (x)) , 



where 



H k (x) 



q ( x k) (x, v)dv + v x k) (x, s) 
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G k (x) 



q x k J(x,v)dv + v x k J(x,s) 



In the numerical experiments the integrals were approximated by a quadrature for- 
mula. It should also be mentioned that although from a theoretical standpoint, the 
parameter s* can be chosen arbitrarily, in the numerical experiments it was chosen to 
provide the best accuracy of reconstruction of some benchmark refraction coefficients. 
The function v(x,s) at s was updated as follows. We first solve numerically the two- 




FlG. 6.5. Reconstructions of the relative electrical permittivity of a low- contrast land mine 
embedded in clutter at x = 0.2 without (asterisks) and with (bullets) refinement. The original 
distribution is shown by a solid line. 

point problem (|6 . T[) - (16.9[) with the (k + l)th successive approximation rik+i{x) of the 
refraction coefficient. Let u^ k+1 ^ (x,s) be such a solution. Then the updated function 
v(x, s) and its two derivatives were computed as 



(x,s) 



(x,s) 



ln(u^ k+1 \x,s)- 


- exp (— sx)) 


X 




s 2 




+ -, 

s 




_(fc + l)/ _N 

U x '(X,S) 


— se~ sx 


1 






- exp (— ~sx)) 


+ =, 

s 




V<xx (^, <s) 


+ s 2 e" SI 


J 


— (k-\-l) f —\ — _ sr 

Ux (x, s) — se 



-1 2 (u( fc+1 ' (x, s) + exp (— sx)) \ s (u( fc+1 ) (x, s) + exp (—sx)) 



Once the inequality 



\\dxU nk (0,S) - <p(s)\\ L2{ s,s) < S 



was fulfilled, the iterative process was terminated. The approximate solution rik Bto (x) 
was refined as follows. We defined the function 

(6.17) u (fcstop) (a;,s0 = u (x, s*) exp [slv ikst " p) (x, s,)] 
and minimized the functional 

(6.18) Jx(n) = ||B B (s,«0-^*^ ) (^«*)lli a (o,i)+^IK«)-"fc« v (a!)lli a (o,i)- 

Since the functional J\(n) is strictly convex, the Powell method with the initial ap- 
proximation n\(x) was used to find its minimizer. 
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Fig. 6.6. Reconstructions of the relative electrical permittivity of the model clutter without 
refinement (bullets). The original distribution is shown by a solid line. 
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Fig. 6.7. Results of recovering the sound speed profile from the data obtained on the see surface: 
reconstructions without (asterisks) and with (bullets) refinement. The original sound speed profile 
is shown by the solid line. 

6.1.1. Numerical results: electromagnetic frequency sounding. We sim- 
ulated electromagnetic frequency sounding of dry layered soil, which is usually per- 
formed by ground penetrating radars for the purpose of land mine detection and 
classification. We assumed that the relative magnetic permeability and electrical per- 
mittivity of the background media are equal to 1, and n(x) — y / e(x)/ea. The realistic 
values of n{x) (see, e.g., |32j ) were used in the numerical experiments. The distribu- 
tion of the relative electrical permittivity of a subsurface land mine was simulated by 
a Gaussian curve. Figure ISTTl shows the typical boundary data <^(s) generated by this 
mine embedded in the homogeneous background (air) with e = 1. In Figure l6~Tl it 
is shown the effect of refinement of E\ (x) via fitting the interior data in comparison 
with refinement via fitting the boundary data. In this experiment, the non-perturbed 
frequency sounding data was used, and the regularization parameter A providing the 
best accuracy was chosen 3.6 • 10" 3 . 

To demonstrate the robustness of the proposed algorithm, we simulated the per- 
turbed data by adding the normally distributed noise to fp(s) 

ip(s) = Tp{s) + -fR, 
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where R is the normally distributed pseudo-random vector with mean zero and stan- 
dard deviation of 1, and 7 is the model standard deviation chosen as 

\m\ 

7 = e M' 

where e = \\f> — p\\/\\ip\\. The perturbed data is shown in Figure [5TT1 Figure P3T2I shows 
the mean relative electrical permittivity obtained from samples containing twenty 
five realization of the pseudo-random vector R. Since robustness of the proposed 
algorithm was also confirmed in all other numerical experiemnts, we limit our next 
demonstrations by the results of recovering the material properties from the noiseless 
data. Since the land mines may be filled with the different explosives and they may 




Fig. 6.8. Results of recovering the sound speed profile from the data obtained on the see floor: 
reconstructions without (asterisks) and with (bullets) refinement. The original sound speed profile 
is shown by the solid line. 

be covered by either metallic or plastic materials, their averaged relative electrical 
permittivity may range from 2 to 25 or more units. Figure 16.31 shows the results of 
reconstruction for the different contrasts. In all such cases, the effect of refinement 
was significant. 

In practice, even a near subsurface, where land mines are usually embedded in, 
is not homogeneous. This may be due to variations of sand and clay porosity, water 
saturation, presence of abris, etc. These inhomogeneities lead to unwanted variations 
in GPR signals that are referred to the clutter. We simulated the model clutter con- 
taining a continuously inhomogeneous layer with the relative electrical permittivity 
ranging from 2.4 to 4.75. Such a variation is typical, for instance, for a mixture of 
sand and clay with the variable porosity and water saturation. The results of recon- 
structions of the model clutter are shown in Figure IB~4l Also, we investigated how the 
presence of clutter affects reconstruction of the relative electrical permittivity of land 
mines. The results of reconstructions are shown in Figures 16.51 and 16.61 Clearly, the 
presence of clutter affects significantly the quality of quantitative imaging. Specifi- 
cally, the higher the contrast the lower the quality, though a near surface land mine 
may still be detected and classified. Furthermore, we observed that refinement by 
fitting the interior data did not lead to a significant improvement if a low contrast 
land mine is embedded in clutter. In case of a high contrast land mine in clutter, 
refinement was not even possible. 
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6.1.2. Results of numerical experiments: acoustic frequency sounding. 

Also, we simulated acoustic frequency sounding of layered marine environments. In 
these cases, n(z) — cq/c(z), where c(z) is the sound speed profile, and cq is the 
reference sound speed. In the numerical experiments, we used the values of parameters 
typical for the shallow waters and water-saturated sediments (see, e.g., [18] )■ To 
simulate the sound speed profile in the water column, we used the simplified Wilson's 
formula [31] 



c(z) = 1449+4.6T(z)-0.055T(z) 2 +0.0003T(z) 3 = (1.39-0.12T(z))(S'(z)-35)+0.017z, 



where T(z) — 20 — (3z is the temperature (Celcius), z is the depth (m), and S(z) — 
34 + tz is the salinity (psu). The sound speed profile in the sea floor sediments 
was simulated by a piecewise-constant average compressional wave speed in the range 
from 1685 to 1710 m/s. Inversion was performed for the dimensionless parameter s 
ranging from 1 to 10 on a grid containing 64 nodes. Figure 16.71 shows the results 
of inversion from the surface of the water column. Although the proposed algorithm 
performed nicely for the water column, it did not allow for obtaining the fine structure 
of sediments. The results of inversion from the sea floor are shown in Figure [oT51 In 
this case, the refinement was resulted in a good reconstructed image, though the 
accuracy of reconstruction was diminished as we advanced into the sediment layer. 
The influence of perturbations in the data on the accuracy of reconstruction was very 
similar to the case of electromagnetic frequency sounding. 

6.2. The Slichter-Langer-Tikhonov (SLT) inverse model of electrical 
prospecting. Recovering the subsurface conductivity in a horizontally stratified 
Earth from the observed surface potential was originated by Slichter [20) and Langer 
[13j . In [2 7) Tikhonov established the uniqueness result for their formulation. In 
[17j the SLT inverse model, as well as the logarithmic derivative of the conductiv- 
ity function and Hankel transform was used to justify the local strict convexity and 
smoothness of the residual functional. In this section we apply the proposed algorithm 
for the SLT inverse problem. 

Assume that the upper-half space is filled with a conductive medium, such that 
a = <j(z) > const. > 0, and there exists a thin layer (— e, 0),e > 0, such that its 
conductivity is constant and equals to cr(0) = 1. The lower half-space z < —s is 
supposed to be filled with a dielectric. Introduce the cylindrical coordinates (r, z) and 
consider a point-like current electrode placed at the point (ro,zo) := (0, Zo),Zq € 
(— e, 0). Because of the cylindrical symmetry, the voltage potential V(r, z) generated 
by this current satisfies the following conditions. 



We formulate the following inverse problem. 

Inverse Problem of Electrical Prospecting. Let the function V(r, z) satis- 
fies the boundary value problem i6.19\) - £6.21\) and it is associated with the conduc- 
tivity o~(z). However, the trace V{r, 0) ofV(r,z) on the surface z — is unknown. 
Instead, its approximation ^(r), such that \\V(r, 0) — ^(r)]! < §, is known. Given 
(^(r), 8, <70) L), find an approximation of o~(z) in the finite layer (0, L). 



(6.19) 
(6.20) 
(6.21) 



r 1 (rV r ) r + a 1 (z)(a(z)V z ) z = -S (r) S(z - z ), r > 0, z > -e, 
a(0)V z {r,-e) = 0, 



| (r.z) | — >-oo 



lim V(r, z) = 0. 
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Following [27] . we first apply the zero-order Hankel (Fourier-Bessel) transform 

/•oo 

V(z,0= / V(r,z)J ((r)rdr, 
Jo 

to the problem (|6.19[) - (|6.21|) and introduce the new variables 

(6.22) t(z) = f'*- l (t)dt, 

x = t/T,s = a(L)^T,n(x) = a(z(xT))/o~(L), where a(t) = o~(z(t)) and 

Then, the problem (j6.19|) - (!6.21[) is reduced to the two-point problem 

(6.23) u xx (x, s) — s 2 n 2 (x)u(x 7 s) = 0, ie (0,1), s e [s,s], 

(6.24) u x (0,s) = -S(x-x ), 

(6.25) u x {l,s) + su(l,s) = 0. 

Here, u(x,s) — V(z(xT),s/a(L)T). We notice the similarity of this problem with 
the problem (|6.4p - (|6.6p modelling frequency sounding of layered media. The only 
difference is the surface boundary conditions (|6.24p . which is due to the presence of 
the current electrode. This observation motivates the following reformulation of the 
inverse problem. 

Let the function u(x, s) be a solution of the problem id. 23]) - [6. 25)) that corresponds 
to n(x), but u(0, s) is unknown. Instead, an approximation <f(s), such that \\u(0, s) — 
< S v , where 5 V > and ip{s) — %jj{s / 'a(L)T) 

/>oo 

ip(r) = / ^(r)J Q (Tr)rdr 
Jo 

is given, find an approximation of n{x) in (0, 1). 

Since the SLT model (|6.23[) - (|6.25[) is similar to the frequency sounding model (|6.4[1 - 
(|6.6p . the specific procedure for reconstructing the conductivity is the same as de- 
scribed in the section [6~T1 However, since the travel time transformation (|6.22j) maps 
a uniform grid {zi}™ =0 onto a non-uniform grid {ti}™ =0 , such that 

U = cr~ 1 (z i )h + U-i, t = 0, t n =T, 

one needs to transform the reconstructed refraction coefficient back to the conductiv- 
ity distribution. Let hi be the values of the reconstructed refraction coefficient on an 
arbitrary uniform grid {xi}" =0 . Then, to obtain the values a of the reconstructed con- 
ductivity on {zi}™ =0 , we interpolated the values hia(L) from the uniform grid {£i}?=o 
onto the non- uniform grid {ti/T}™ =0 corresponding to the uniform grid {zi}™ =0 . In 
this case, the "optimal" dimensionless range of the parameter s, i.e., [s,s], was chosen 
as [0.01,0.1]. 

In case of the SLT model of electrical prospecting the refinement procedure did 
not lead to the significant improvement of cr/c stop (z). Therefore, Figures IB31 and [6.101 
show the results of reconstructions without refinement. 
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Fig. 6.9. Inversion of the surface voltage potential (bullets). The original distribution of 
conductivity for the two-layer medium is shown by the solid line. 
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Fig. 6.10. Inversion of the surface voltage potential (bullets). The original distribution of 
conductivity for the four-layer medium is shown by the solid line. 



7. Conclusions. We developed a globally convergent algorithm for quantitative 
imaging in frequency sounding and electrical prospecting of layered media. The con- 
vergence of the algorithm was justified. Compared to the existing globally convergent 
algorithms, the proposed one possesses two new features. These are the iterative 
procedure for the numerical solution of an overdetrmined boundary value problem 
for the second-order nonlinear differential equation and the procedure of refinement 
of an approximate solution. It was demonstrated in the numerical experiments that 
these novelties may improve, sometimes significantly, both the spatial and contrast 
resolution of imaging. 
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